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Abstract 

The evolution of a pile of granular material is investigated by molecular 
dynamics using a new model including nonsphericity of the particles 
instead of introducing static friction terms. The angle of repose of 
the piles as well as the avalanche statistics gathered by the simulation 
agree with experimental results. The angle of repose of the pile is 
determined by the shape of the grains. Our results are compared with 
simulations using spherical grains and static friction. 

PACS. numbers: 05.60. +w, 49.29.-i, 46.10,+z, 81.35.+k 
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1 Introduction 



There are many interesting phenomena observed in the dynamics of granular 
materials which have been studied by many authors and by various methods 
(e.g. 0J0] an d references therein). One of the most interesting phenomena is 
the evolution of a pile of granular material. When sand grains are dropped 
on the top of a pile its slope does not depend on the number of grains as 
long as the heap is not too small. The avalanches going down the surface of 
the heap have no characteristic size, they are only limited by the size of the 
entire pile. The size of the avalanches as well as the time interval between 
successive avalanches fluctuate irregularly. Their spectra decay possibly due 
to a power law On large scales, however, there was measured a differ- 

ent behaviour ||. Recently Rosendahl et. al. found that the predominant 
number of avalanches is power-law-distributed but that these avalanches are 
supplemented by periodic significantly larger avalanches. Dependent on the 
material the slope of the surface of a heap can vary, hence one observes waves 
creeping up the pile JTD[ . 

There are different numerical methods to simulate the static and dy- 
namic behaviour of macroscopic amounts of granular media. Cundall and 



Strack JTTJ] introduced a model for the interaction of grains which is widely 
used in molecular dynamics simulations. Other simulations base on cellular 
automata [12| , the results led to the idea of self organized criticality . 
Other numerical investigations base on the Boltzmann lattice gas ||14|| , on 



Monte Carlo simulations (e.g. (T^][|T^]) or on random walk models [[17 



In this paper we apply a new model for the simulation of the evolution 
of sandpiles using molecular dynamics as introduced in Many of the 

numerical investigations base on molecular dynamics (e.g. p |]lTl||T9| -p6|), 
however, most of them use spherical grains. To simulate static friction be- 
tween the grains one introduces a force which is due to the Coulomb law, 
i.e. the spheres slide on each other for the case that the shear force between 



the grains Fs is larger than /i • Fn where Fn is the force in normal direction, 
H is the Coulomb friction coefficient. This friction force is derived from the 
phenomenological idea about friction between macroscopic bodies. Its mi- 
croscopic origin is still hardly understood f27j. In our molecular dynamics 
simulations we use nonspherical particles instead of spheres. The particles we 
use here are elastic but not stiff, they deform during collisions and dissipate 
energy. 

Recently Gallas and Sokolowski suggested a model for nonspherical grains 
where two spheres were connected by a stiff bar |28]. Their results agree 
qualitatively with ours. 

As we will demonstrate below the results of the simulations using non- 
spherical grains agree better with experimental investigations than tradi- 
tional models. 



2 The Model 

In our simulation each of the grains k consists of five spheres, four of them 
with radii rf^ are located at the corners of a square of size The fifth 

sphere is situated in the middle of the square, its radius is chosen to touch 
the others 

Each of the spheres is connected to all of its neighbouring spheres which 
belong to the same grain by damped springs (fig. [1]). 

Figure 1: Each of the 'particles k consists of four spheres with equal radii 
at the corners of a square of size and one sphere in the centre. Its radius 
is chosen to touch the others. There are forces due to a damped spring acting 
between the spheres the grain consists of. 



Between each two spheres i and j which might belong to the same particle 



or to different particles acts the force F™: 



1-U = { ' ' ' A (2) 
otherwise 

with the normal force 

F£ = Y ■ (r { + Tj - \Si - Xj\) - 7 • m\y ■ \x { - Xj\ (3) 
and the effective mass 

eff rrii ■ rrij 

™4 = ■ 4 

J rrii + rrij 

Y and 7 are the Young modulus characterizing the elastic restoration of the 
spheres and the phenomenological damping coefficient. The terms Xi, X{ and 
rrii denote the current position, velocity and mass of the i-th sphere. The 
normal force consists of an elastic restoration force which corresponds to 
the microscopic assumption that the particles can slightly deform each other 
and a term describing the energy dissipation of the system due to collisions 
between particles according to normal friction and plasticity. 

Moreover each pair of neighbored spheres i and j which both belong to 
the same grain k feel the force 



F* 



s P 



a- (C\f - Xj\) -js P - ^ • Wi - xj\ ~ — (5) 
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with 



' ^ else . 

This force is due to in internal damping with spring constant a and damping 
coefficient 7s p acting between particles of the same grain. 

To compare our results using the nonspherical grains defined above we 
describe now how to introduce static friction in the simulation using spheres 



if i, j lie at the same edge 



as it was done first by Cundall and Strack [D]] and modified by Haff and 
Werner j0|. Here each pair of particles i and j with radii and 77 feel the 



force 



p = 1 V \x l -S J \ V \1 ; |x*i-5j| - ^ 

otherwise 



with the force in normal direction 



F£ = Y ■ (r { + rj - \Si -Xj\) + 7 • mly ■ |fj - Xj\ (8) 

and the shear force 

F§ = min{- 7s • m$f ■ \^f\ ,ii-\Fg\] • (9) 

The relative velocity of the particle surfaces results from the relative velocity 
of the particles and their angular velocities Qi and Qj 

v[f= (Si-'x^+n-ili + Tj-tlj . (10) 

The parameters 7 and 75 stand for the normal and shear friction coefficients, 
/i is the Coulomb friction parameter. Eq. (Q) takes into account that two 
grains slide upon each other when he shear force between them exceeds a 
certain value, otherwise they roll (Coulomb relation). 

The sizes of the nonspherical particles as well as the radii in the 
case of spherical grains were Gauss distributed with mean value L( fc ) . 

3 Results 

In our simulation we investigate the evolution of a sand pile by consecutively 
dropping nonspherical particles on the top of the heap. The surface of the 
platform upon which the heap is built up consists of spheres with radii r- p ^ 
with mean value = to simulate a rough surface. After each drop- 
ping a particle on the heap we waited until the velocities of all particles of 
the entire heap faded away, i.e. until they were smaller than a very small 
constant. 

For the parameters we used in the simulations: Y = 10 A kg/s 2 , 7 = 
1.5 • 10 4 s-\ a = 10 4 kg /s 2 , 7 5p = 3 • 10 4 s~\ /rf ] = 4, W) = 3 mm. 
The results were compared with simulations using spherical grains, here we 
assumed the same values and 75 = 3 • 10 4 fi = 0.5 and r\ — 3 mm. 



For the calculation of the evolution of the particle positions we applied 



a Gear predictor corrector algorithm [2G] of sixth order. For the case of the 
simulation using spherical grains we have further to calculate the rotation 
of each sphere. Since the acceleration of the particle surface due to particle 
rotation is much smaller than due to particle movement we applied therefore 
a Gear algorithm of fourth order only. 



3.1 Angle of Repose 

When building up a heap experimentally one finds that the slope $ (angle 
of repose) does not depend on the particle number. It depends on the raw 
material and lies typically between 20° and 30°. Bretz et. al. found $ ~ 
25° §. 

Figure |] shows snapshots of simulated piles of different sizes. The size 
of the piles have been scaled so that the widths of all piles appear equal in 
the figures. The first two piles consist of N = 300 (a) and N = 1100 (b) 
nonspherical particles. The slope is approximately the same for both figures, 
i.e. it does not depend on the number of particles. Figs, (c) and (d) show 
piles of N = 400 and N = 1800 spherical particles. In the case of spherical 
grains the heap dissolves under gravity, the slope depends on the number of 
particles. 

Figure 2: Snapshots of simulated piles. Heaps (a) and (b) consist of N = 300 
and N = 1100 nonspherical grains, heaps (c) and (d) of N = 400 and N = 
1800 spherical particles. In the latter case the slope depends on the number 
of particles. 



In the simulations using nonspherical grains we measured $ « 21°, this 
value is in agreement with experimentally found data ||. Fig. |3| shows the 
evolution of the slope $ of the pile during one run over the particle number 
for spherical and nonspherical grains. For nonspherical grains the slope $ 
fluctuates due to avalanches of different size. Except for very small heaps the 
average slope does not depend on the particle number in accordance with the 



experiment. For the case of spherical particles the slope decays with rising 
particle number. 

Figure 3: The evolution of the slope of the pile during one run over the 
particle number for both spherical and nonspherical grains. For nonspheri- 
cal grains the mean angle does not depend on the particle number while it 
decreases for spherical particles. 



Since the surface of the pile is not a smooth line we need a special pro- 
cedure to calculate the slope $. Fig. |] shows a pile with an arbitrarily 
shaped heap. The shape of the pile is described by the function h(x) with 
x G [0, x max ], x cm and y cm denote the coordinates of the centre of mass of 
the heap. The height H of an ideal triangle which has the same baselength 
and volume like the heap would be 

H = f h(x)dx . (11) 

Because the centre of mass of a triangle is the intersection point of all three 
median lines we calculate the baselength B of the triangle using the intercept 
theorem 

2 ■ H _ H - y cm 
B ~ x ' 1 ' 

Hence we find for the slope 



$ = arctan ( H Vcm ) . (13) 



2 ■ x cm 

Eq. ( |T3"D gives a good approximation for the slope provided that the shape 
of the heap is close to a triangle. The snapshots in fig. [2] show that this 
precondition is given in our case. 



Notice that there are two possibilities of defining the angle fl3(J , the angle 
of repose and the angle of marginal stability. For simplicity we used in this 
work the time average of the angle observed during the simulation. 



3.2 Avalanches and Mass Fluctuations 



Theoretical as well as experimental investigations PJ-|@][r3| led to the hy- 
pothesis that the size of the avalanches, i.e. the mass fluctuations of the pile 



Figure 4: Schematic plot of the procedure used to determine the slope of a 
pile which surface is not a smooth plane. 

as well as the distribution of the time intervals between successive avalanches 
of sandpiles can be described by the self organized criticality-model. In order 
to investigate the size of the avalanches and the distribution of the time be- 
tween each two successive avalanches we modified the setup of the simulation 
and monitor the fluctuation of the mass of a heap of definite size. Instead 
of an infinite platform we use a platform of a finite length P above which the 
heap is built up. In fig. |5| is drawn a part of the time series of the mass mh 
for fixed P = 820 mm w 273 L^ k \ The mass fluctuates irregularly due to 
avalanches of different size going down the surface of the heap. Since after 
dropping a grain on the heap we wait until the motion of all grains fades away 
the natural unit for measuring the time is the number of dropping events, 
but not the number of integration steps used in the simulation. Since we 
want to make statistics for all avalanches but not only for those which cause 
a change of the mass mh of the heap we monitor the change of the slope of 
the heap, i.e. the change of the centre of mass, during an avalanche instead 
of the change of the mass. Otherwise we would cut off the lower part of the 
spectrum, most of the small avalanches do not reach the bottom of the pile. 
Fig. | shows the time series of the avalanches. Fig. ||a displays the changes of 
the slope, fig. |B|b the changes of the mass. Both figures resemble each other 
but in the latter case many of the small avalanches have been cut off. This is 
due to small avalanches which do not reach the bottom and therefore do not 
change the mass of the heap. The log-log-plot of the spectrum (fig. [7]) shows 
that the size distribution of the avalanches might follow a power law. For 
the exponent H(Na) ~ {Na} 171 we measured m « —1.4. In the experiments 
was found m ~ —2.5 0] and m « —2.134 || respectively. The spectrum of 

Figure 5: Total mass of a pile of nonspherical grains on a platform of finite 
length P = 820 mm ~ 273 L( k h The mass fluctuates irregularly due to 
avalanches of different size. The time is measured in dropping events. 



Figure 6: Time series of the changes of the slope (a) and of the total mass 
(b) of the heap due to avalanches. Many of the small avalanches do not cause 
a change of the mass of the heap, since they stop before reaching the bottom. 

Figure 7: Size distribution of the avalanches. The line displays the function 
h(N A ) ~ (N A )~ 1A . 

the time intervals between each two successive avalanches was experimentally 
investigated too. Fig. |8| shows the spectrum of the distribution of the time 
intervals between two successive avalanches which we found in the simula- 
tion. We are not convinced that the decay of the spectrum in our simulation 
follows a power law too. 

Figure 8: Spectrum of the time intervals between successive avalanches. 



3.3 Influence of the Particle Sphericity 

Dependent on the ratio L^/r^ the grain k is more similar to a square or 
to a sphere. Hence we may investigate the influence of the particle shape on 
the angle of repose of a pile. We define the shape S of a grain 



>cc 



DC 

S = 1 _ fhBEL (14) 



fee 
L max 



where R^ in and R^ ax are the extremal values of the distance between the 
convex cover of the nonspherical grain and its central point (fig. |9|). In terms 
of /rf } we write 

/ i L (k) i L (k) \ 

1 

Rmax = ~^J^ (k) 1 (16) 

Fig. [10] displays the shape value S as a function of /r\ . The function 
S reaches its maximum S m for a particle which convex cover is most similar 
to a square ((L^ /r^)s m = 9.66). For each S < S m there are two different 
particle shapes with the same shape value. Since L^/r^ > 2 there are no 



Figure 9: The shape S is determined by the extreme sizes of the convex cover 
of the grains R^ in and BZ 



max ' 



particles with S < 0.172 for lower L^/rj ratio. For the limit L^/r| — > oo 
the grains are spheres and S becomes zero. 

Figure 10: The shape value S as a function of iP^/rf^. Each value S with 
S < S m corresponds to two different particle shapes with the same shape 
value, one with lower and one with larger L^/rj ratio. 

To investigate the influence of the shape S of the grains on the result 
of the simulations we have to scale the density of the material the grains 
consist of to ensure that the total mass of each grain remains constant 
independent on its shape (po - ► P when So — > S, L^ k > = const., = const.) 

(2-S) 2 45 2 -45 +2 , L« ^ ( L« 



P° ' (2-S ) 2 4S 2 -4S+2 iU1 r W - I r « 

p(S) = { , , 1 ' m (17) 

1 4S 2 165 , ( 2 + (l2y / 2-24)5o+15-10v / 2 L (k) ( L (k)\ 

P ° ' 45j 165 2 + (l2^-24)S+15-10 v / 2 tOT 7^ V^F/ 5m ' 

The influence of the shape of the grains to the slope of the heap is drawn in 
fig. [TT| As mentioned above each S corresponds to two arguments L^>/r\ . 
Therefore we denote the slope for grains with L^/r} < (IS® Jr % )s m by 
© and by + for L^/rf } > (L^ /rf ] ) Sm . As expected the slope $ reaches 
its maximum for particles with maximum shape value S = S m since they 
are most similar to squares. The dashed line marks the inclination $ sp we 
found for a heap of spheres, which corresponds to S —>■ 0. The slope values 
we measured simulating nonspherical grains lie between the slope of a pile of 
spheres and a pile of grains with S = S m . 

Figure 11: Slope $ of a heap over the shape value S for grains with L^/r^ < 
{L^/r\ k) ) Sm (0) andL^jrf > (L^/r\ k) ) Sm (+). The dotted line leads 
the eye to the function $ = 130 ■ S + const. The dashed line displays the 
inclination observed in simulations with spherical particles. 



The calculation for the data shown in fig. [IT] are extremely time consum- 
ing. For this reason we are not able to present more data. 



4 Conclusion 



Using molecular dynamics simulations we investigated the evolution of a sand 
pile and the statistics of the avalanches going down. Simulated piles using 
nonspherical grains remain stable under gravity, independent on their size the 
slope is approximately the same. In the simulation we found for the average 
slope $ fa 13° . . . 23° dependent on the shape of the grains. In experiments 
was found $ w 25° || . For the case of simulations using spheres and friction 
forces due to the Coulomb law the slop of the heap decreases with increasing 
particle number. 

During the simulation the mass of the heap which is built up upon a finite 
platform varies irregularly. For the statistics of the avalanches we used the 
fluctuations of the slope of the heap, using the fluctuations of the mass leads 
to wrong results since most of the small avalanches do not reach the bottom 
of the heap and therefore do not change its mass. In our simulations the 
avalanches are power-law-distributed, for the exponent we found h(N^) ~ 

(n a )- 1a . 

Our model allows for changing the shape of the grains. We investigated 
the influence of their shape to the slope of the heap of constant size. We 
found that the slope reaches its maximum for grains which convex cover is 
most similar to a square. The minimum slope appears for the case that the 
particles are spheres. The slope of piles built of other nonspherical particles 
lies in between the minimum and the maximum. 
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